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Abstract Evidence for an anomalous annual periodicity in certain nuclear decay 
data has led to speculation concerning a possible solar influence on nuclear 
processes. We have recently analyzed data concerning the decay rates of 36 C1 
and 32 Si, acquired at the Brookhaven National Laboratory (BNL), to search for 
evidence that might be indicative of a process involving solar rotation. Smoothing 
of the power spectrum by wcightcd-running-mean analysis leads to a significant 
peak at frequency 11.18 yr -1 , which is lower than the equatorial synodic rotation 
rates of the convection and radiative zones. This article concerns measurements 
of the decay rates of 226 Ra acquired at the Physikalisch-Technische Bunde- 
sanstalt (PTB) in Germany. We find that a similar (but not identical) analysis 
yields a significant peak in the PTB dataset at frequency 11.21 year -1 , and a 
peak in the BNL dataset at 11.25 yr -1 . The change in the BNL result is not 
significant since the uncertainties in the BNL and PTB analyses are estimated 
to be 0.13 year -1 and 0.07 year -1 , respectively. Combining the two running 
means by forming the joint power statistic leads to a highly significant peak at 
frequency 11.23 yr -1 . We comment briefly on the possible implications of these 
results for solar physics and for particle physics. 
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1. Introduction 



For the last one hundred years, scientists have believed that the decay rate of 
each radioactive element is constant, unaffected by any environmental process. In 
1930, Rutherford and his colleagues wrote "The rate of transformation of an ele- 
ment has been found to be a constant under all conditions" ( Rutherford et al., 1930[ )[ 
. In the latter half of the 20th century, however, there were indications that it is 
possible to alter the decay rates of certain isotopes ( |Emery, 1972[ |Hahn, 1976| 
|Dostal et al, 1977] ) by physical or chemical processes. 

We have recently found strong evidence ( Jenkins et al., 2009b Fischbach et al, 2009| | 
|Javorsek II et al., 2010[ ) for an annual periodicity in decay data acquired at 
the Brookhaven National Laboratory [BNL; ( |Alburger et al., 1986 )] and at the 
Physicalisch-Technische Bundcsandstalt [PTB; ( jSiegert et al., 19"98| ] [although 
there appears to be a difference between early PTB data and more recent data 



(Schradcr, 1992 Schrader, 20101]. The BNL experiment monitored the decay 
rates of a2 Si and^Cl over the interval 1982 to 1989. (Alburger et al. preferred 
to analyze only data from 1982-1986, but we have preferred to analyze the 
entire data set.) The PTB experiment monitored the decay rates of 152 Eu, 
154 Eu, and 226 Ra from 1983 to 1998. These analyses, and an apparent change 
in the measured decay rate of 54 Mn during a solar flare in December 2006 
( | Jenkins and Fischbach, 2009a[ ), have led to the suggestion that some nuclear 
decay rates are influenced by particles (possibly neutrinos) or fields emanating 
from the Sun fl Jenkins et al., 2009b] |Fischbach et al., 2009| . 

As one would expect, the suggestion that nuclear decay rates may be variable 
has not gone unchallenged: 



i) 



ii ) 



in) 



Scmkow et al. (2009) and others have suggested that these fluctuations have 
their origin in environmental or systematic processes. However, detailed in- 
vestigation shows that the results of our BNL and PTB analyses cannot be ex- 
plained by variations of temperature, pressure, humidity, etc. ( | Jenkins et al., 2010[ ).[ 
Norman et al. (2009), have re-examined data from several nuclear-decay 
experiments, finding no evidence for a correlation with Sun-Earth distance. 
However, our collaboration has re-analyzed the Norman group's data, which 
Norman and his collaborators generously provided, and we have detected an 
annual periodicity, small in amplitude but with the same phase as we have 
found in the BNL and PTB datasets. 

Cooper (2009) has analyzed data from the power output of the radioisotope 
thermoelectric generators aboard the Cassini spacecraft, finding no significant 
deviations from exponential decay, but we find no conflict between Cooper's 



results and our results (Jenkins et al, 20101 



To pursue the question of a possible association between solar processes 
and variations in nuclear decay rates, we have recently examined the BNL 
dataset, searching for a periodicity that may be associated with solar rotation 



(Sturrock et al., 2010a|-in particular, with rotation rates that show up in low- 



energy solar neutrino data. We have found that low-energy solar neutrinos, 



as detected by the Homestake ( |Davis et al, 1968] |Cleveland et al, 1998Q and 
GALLEX ( |Anselmann et al, 1993[|Anselmann et a~1 995; Hampcl ~et al, 1996[ | 
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Hampel et al, 1999[ ) experiments, exhibit a periodicity at about 12 year 1 ( Sturrock, 2008[ [ 



Sturrock, 2009), which we find also in ACRIM irradiance data ( |Willson, 1979[ 



Willson, 2001 ). Since this frequency is significantly lower than the dominant 



synodic rotation frequency of the radiative zone, which is about 13 year -1 (14 
year -1 sidereal; ( Schou et al, 1998 1), we have suggested that it may be related 
to the synodic rotation frequency of the solar core (the rotation rate of the core 
is ill-determined at this time). 

In our recent analysis of BNL data, we find a strong peak at 11.18 year -1 
dSturrock et al, 2010a[ ). The purpose of this article is to carry out a similar anal- 
ysis of PTB data to see whether these data also show evidence of a periodicity 
at a frequency lower than the radiative-zone rotation rate. 

The preparation of the data and the power spectrum analysis arc presented 
in Section 2, where we focus on a "search band" of 10-15 year -1 . This band 
extends high enough to cover the equatorial synodic frequencies of the radiative 
and convection zones ( |Schou et al, 1998[ ), and extends well below the frequency 
of the strong peak in the BNL power spectrum. In order to assess the probable 
errors of our frequency estimates, we note that the Fourier transform of a pure 
sinusoidal signal of length T at frequency v D peaks at v D and falls to zero at 
v Q ± 8v, where 5v = 1/T. Since the durations of the BNL and PTB datasets are 
7.80 years and 14.84 years, we adopt the generous error estimates 0.13 year -1 
and 0.07 year -1 , respectively. With this guideline, we find a feature in the PTB 
power spectrum at frequency 11.29 ± 0.07 year -1 , close to the frequency of the 
major peak in the BNL power spectrum. We examine the significance of this 
periodicity in Section [31 using both the standard shuffle test and a variant of 
this test that we call the "shake" test. 

Since the power spectrum is more complex than we would expect of a stable, 
rigid rotator, we adopt in Section [3] the procedure we used in our analysis of 
BNL data-that of forming the weighted running mean of the power spectrum, 
which leads to a prominent feature at frequency 11.25±0.07 year -1 . In Section 
21 we examine the correlation between the BNL and PTB datasets by forming 
the "joint power statistic" derived from the BNL and PTB weighted-running- 
mean power spectra. This shows one dominant feature at frequency 11.23±0.07 
year -1 . We evaluate the significance of this feature by again using the shuffle 
test. We discuss these results in Section [5j 



2. Power Spectrum Analysis 

The PTB dataset comprises 1966 measurements of the decay rate of 226 Ra in the 
time interval 1983.855 to 1998.809, where dates are measured in what we call 
"neutrino years," which have proved useful in analyzing neutrino data. Dates 
in "neutrino days" are counted from January 1, 1970, as day 1, and dates 
in "neutrino years" are given by 1970 + (neutrino days)/365.2564. We have 
prepared the data for time-series analysis by dividing the decay count rates 
by the counts expected on the basis of a purely exponential decay, using the 
mean decay rate determined from the best fit to the data. We also remove a 
few "outliers," namely datapoints for which the measurements differ from the 
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mean by more than three sigmas. The resulting normalized count rate is shown 
in Figure [TJ in which the annual variation is obvious. 




1&S4 I9£6 1968 1990 1992 1994 1996 1996 2000 
Date 



Figure 1. PTB decay-rate measurements, corrected for mean decay rate and normalized to 
mean value unity. 



We have next performed a power-spectrum analysis of the data shown in Fig- 

1 using a likelihood procedure ( |Sturrock et al, 2005a[ )| 



ure[TJ in frequency steps of 0.01 year 
that is equivalent to the Lomb-Scargle procedure ( |Lomb, 1 976 Scargle, 19821 
The result is shown in Figure [2] There is, as expected, a huge peak at 1.00 
year -1 , with power S = 507, due to the annual modulation which is obvious in 
Figured] Between 0.01 year -1 and 2.11 year -1 , there arc 13 peaks with power 
20 or more. We also note that there are strong peaks near 51 and 53 year -1 (not 
shown in this figure), which are obviously due to aliasing of the peak at 1.00 
year -1 due to a weekly pattern in data acquisition. 

Our main goal in this article is to search for features in the power spectrum 
that might be related to solar rotation. We shall adopt a search band of 10-15 
year -1 that is wide enough to include periodicities found in previous analyses of 
related data-neutrino data and irradiance data ( |Sturrock~ 2008; St urrock, 2009| ), 
and BNL decay data ( Sturrock et al., 2010a ). It is also wide enough to include 
the band of synodic rotation rates of an equatorial section of the convection 
and radiative zones ( |Schou et al., 19 98). With this goal in mind, it is prudent 
to remove the influence of the very strong low-frequency modulations, which 
can cause aliasing and otherwise adversely influence the power spectrum. We 
find we can suppress most of the low-frequency modulation without appreciably 
degrading modulation within or above the rotational band by the following 
procedure. 

We form 11-point running means of the normalized count rates, and then sub- 
tract the running- mean values from the normalized rates. The resulting modified 
count rates are shown in Figure [3] A likelihood analysis of this time series yields 
the power spectrum shown in Figure 0J The principal feature is at frequency 
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Figure 2. Power spectrum of PTB data, formed by a likelihood procedure. 



11.29 ± 0.007 year 1 with power S = 15.20, but there is a cluster of peaks in 
that vicinity. 




1992 WW 1SS3 158$ 1990 1992 S994 1996 199a 2000 
Data 

Figure 3. PTB decay-rate modified count rates measurements. 



Noting there are 48 peaks in the band 10-15 year we find that the standard 
formula for the false-alarm probability flScargle ,~1 982) 



FAP = 1 - ( 1 - e" 



(1) 
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Figure 4. Power spectrum of modified PTB data, formed by a likelihood procedure. The 
principal peak is found at 11.29 year -1 with power S = 15.20. 



yields the estimate 1.2x 10~ 5 . However, as we have recently pointed out ( jSturrock and Scargle, 2010bp ,[ 
if one is over-sampling (as in the present case) , one needs to replace S in Equation 
[T]by S — 1. With this change, Equation [1] yields the revised estimate 3 x 10~ 5 . 

In order to obtain a more robust significance estimate, we have carried out 
the familiar shuffle test ( Bahcall and Press, 1991 1. We find that, of 10,000 such 
simulations, none has a power as large as 15.20 in the search band. The distri- 
bution is shown as a histogram in Figure [5] and in logarithmic form in Figure [51 
A projection of the distribution shown in Figure [5J indicates that the probability 
of obtaining by chance a power as large as or larger than 15.20 is only 5 x 10~ 5 , 
close to the FAP value. 

We have also applied the shake test ( |Sturrock et al, 2010a| ) that involves 
small random time displacements. The results are very similar to those obtained 
from the shuffle test. 



3. Frequency- Averaged Power Spectra 

As in our analysis of BNL data, we see from the power spectrum shown in 
Figure 4 that there is actually a cluster of peaks in the vicinity of 11 year -1 , 
indicating that the time series should not be attributed to a single stationary 
periodic process. The combination of periodicities may be due to the influence of 
different sources with different characteristic time scales, for example, different 
rotation rates. For this reason, we adopt a procedure for evaluating a complex of 
peaks in the search band, rather than focusing on only one peak. This procedure 
is that of forming a weighted running mean of the power spectrum. 

We denote by Vj the sequence of frequencies (with spacing 0.01 year -1 ) at 
which the power is evaluated, and by Sj the power sequence, where j = 1, . . . , N v . 
We then form the sequence 
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Figure 5. Histogram display of the result of 10,000 shuffle simulations of the modified PTB 
dataset. None of the simulations yields as big a power (15.20) in the search band (10-15 year -1 ) 
as the actual dataset. 
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Figure 6. Logarithmic display of the result of 10,000 shuffle simulations of the modified PTB 
dataset. A projection of the curve indicates that there is a probability of obtaining by chance as 
big a power (15.20) in the search band (10-15 year -1 ) as the actual dataset is about 5 X 10 — 5 . 



S = J^W(j-k)S k . (2) 

k 

We have adopted 
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in which C is chosen so that the sum of the weighting terms W is unity: 



(3) 



C 



> cos — — 

— in v 



(4) 



We have adopted m = 50, which is equivalent to forming running means of the 
power over intervals of width 1 year -1 . However, we find that 



1/2 



W{j)f 



0.435to, 



(5) 



so that, for m = 50, the rms frequency deviation is only 0.22 year -1 . 

We show in Figure [7] the result of applying this smoothing operation to the 
power spectrum shown in Figure H]. We find that the smoothed power spectrum 
has its principal feature at 11.21±0.07 year -1 with peak power 4.51. 




10 

Frequency (yeei~ 1 ) 

Figure 7. Plot of the 101-point weighted running means formed from the power spectrum 
shown in Figure[4] The peak occurs at 11.21 year -1 with weighted-running-mean power 4.51. 



Since the running means are not distributed exponentially, it is necessary 
to carry out the shuffle test or a similar test to obtain a significance estimate. 
We again use the shuffle test, re-computing the power spectrum and then the 
running mean of the power spectrum many times, keeping the original times and 
the original count rates, but randomly sorting one of these datasets. We have 
carried out this procedure 10,000 times, determining the maximum weighted 
running mean of the power in the rotational search band for each simulation. 
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The histogram formed from these maxima is shown in Figure |S] We see that 
none of the random simulations yields as large a value of the running mean of 
the power in the search band as the actual value (4.51). 

2000 i 1 1 1 1 , 

1800 
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Figure 8. Histogram of the wcightcd-running-mcan powers computed for 10,000 shuffle 
simulations of the modified PTB data. None is as large as the actual value (4.51). 

These results arc shown in a logarithmic display in Figure [jj, which also shows 
an extension of an empirical curve derived from a projection of the last 1,000 
points. This projection leads us to a p-value of 10 -7 ' 7 , implying that we could 
expect to obtain a peak as large or larger than the actual peak only once in 
about 2 x 10 s simulations. 

We have also examined the running-mean power spectrum by means of the 
shake test, again carrying out 10,000 simulations. We find that the results are 
very similar to those obtained by the shuffle test. 



4. Comparison and Combination of BNL and PTB Data 

We now derive weighted-running-means of the BNL power spectrum. To be 
consistent with our analysis of PTB, we here again adopt m = 50. We show the 
result, together with the PTB running- mean power spectrum for comparison, 
in Figure [TU1 We find that the BNL curve has a peak at 11.25 year -1 with 
uncertainty 0.13 year -1 and running-mean power 7.74. This is indistinguishable 
from the peak in the PTB curve at 11.21 year -1 that has uncertainty 0.07 year -1 . 
Taking into account the uncertainty of 0.13 year -1 , it is also indistinguishable 
from our earlier estimate of 11.18 year -1 ( jSturrock et al., 2010ap . 

It is useful to examine the degree of correlation between the BNL and PTB 
power spectra. A convenient procedure is to form the "joint power statistic" 
( jSturrock et al., 2005b ). For the present situation, in which we are considering 



only two power spectra, we first form the geometric mean: 
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Figure 9. Logarithmic plot of the weightcd-running-mcan powers computed for 10,000 shuffle 
simulations of the modified PTB data. A projection indicates that there is a probability of 
only about 2 X 10 — 8 of obtaining by chance a values as large as the actual value (4.51). 
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Figure 10. Plots of the 101-point weighted running means formed from the PTB power 
spectrum (blue) and from the BNL power spectrum (magenta). The two peaks are found at 
11.21 year -1 and at 11.25 year -1 , respectively. 



X = (S 1 S 2 ) 1/2 . (6) 
The joint power statistic is then defined by 



J 2 = - In (2XKi (2X)) (7) 
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where K\ is the Bessel function of the second kind. This statistic has the useful 
property that if each power is distributed exponentially, so that 



P(S)dS = e- b dS, (8) 

then J satisfies the same distribution. We find that the following simple expres- 
sion offers a close approximation to that of Equation [7J 



1.943X 2 

J2A ~ 0.650 + X- (9) 

We now compute the joint power statistic for the running-means that we have 
calculated for the BNL and PTB power spectra. This offers a useful comparison 
with the running-mean of each power spectrum. We show the result of this 
calculation in Figure 11. We find a prominent feature at frequency 11.23 ± 0.07 
year -1 with J = 10.34. 




°0 5 10 15 so 

Fiequensy (yew -1 ) 

Figure 11. The joint power spectrum formed from the weighted running-mean BNL and PTB 
power spectra. The peak is found at 11.23 year" 1 with J =10.34. 

We again use the shuffle test , computing 10,000 simulations of the calculations 
that led to Figure QTJ For each simulation, we retain the actual measurements 
and the actual times for each dataset, but re-assign them randomly. For each 
simulation, we compute the two power spectra, form the running means, and 
then find the maximum value of the joint power statistic in the search band 
10-15 year" 1 . 

We show the distribution of these maximum values in histogram form in 
Figure [T^] , and in logarithmic form in Figure 1131 The maximum value found in 
these 10,000 simulations is 3.67, whereas the value obtained from actual data 
is 10.34. Figure Q~3] shows an extension of an empirical curve derived from a 
projection of the last 1,000 points. This projection leads us to a p- value of 10 -172 , 
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implying that we could expect to obtain a peak as large or larger than the actual 
peak only once in about 10 simulations. 
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Figure 12. Histogram of the weighted-running-mean powers computed for 10,000 shuffle 
simulations of the joint power statistic formed from the BNL and PTB running-mean power 
spectra. None is as large as the actual value (10.34). 
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Figure 13. Logarithmic plot of the joint power statistic computed for 10,000 shuffle simula- 
tions of the joint power statistic formed from the BNL and PTB weighted-running- mean power 
spectra. A projection indicates that there is a probability of only about 10~ 17 of obtaining by 
chance a value as large as the actual value (10.34). 
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5. Discussion 



One must be cautious about accepting a p- value as small as 10 -172 at face value. 
Nevertheless, there is clearly strong evidence that both BNL and PTB decay- 
rate measurements arc influenced by a common periodicity with a frequency of 
about 11.25 year -1 . This is too low a frequency to reflect a rotational modulation 
that originates in the radiative zone or a near-equatorial region of the convection 
zone. This leaves only the innermost region of the Sun, with normalized radius 
of approximately 0.2 ( |Garcia et ai, 2007] ), for which we do not yet have reliable 
rotation measurements. 

However, the possibility that decay-rate variations are due to solar neutrinos 
provides an independent reason for focusing attention on the core, since that is 
the location of nuclear reactions that produce neutrinos. It is also relevant to 
recall that we have previously found evidence that low-energy solar neutrinos 
exhibit a periodicity at or near 11.85 year -1 ( Sturrock, 2008 Sturrock, 2009 1, 
lower than the estimated synodic rotation rate of the radiative zone. These 
considerations suggest that the sub-radiative-zone region of the Sun is not in 
rigid rotation. If this is so, we must expect that neutrinos of different energies will 
have different periodicities. As a consequence, it is possible that the decay-rates 
of different nuclei will be different. 

A possible scenario is that the central part of the core has a synodic rotation 
rate of about 11.25 year -1 , corresponding to a sidereal rotation rate of 12.25 
year -1 . There must then be a radial gradient of rotation rate between the core 
and the radiative zone, that is believed to have a sidereal rotation rate of about 
13.60 year -1 (Sch ou et ai, 1998) ). Such an intermediate region would constitute 
an "inner tachocline," which might lead to dynamo action, similar to the dynamo 
action attributed to the outer tachocline ( |Krause et ai, 1992[ ). The periodicity 
of 11.85 year -1 , detected in low-energy-neutrino data and total solar irradiance 



data (Sturrock, 2008 Sturrock, 2009), corresponding to a sidereal rotation rate 



of 12.85 year , would originate somewhere in the inner tachocline. 

A rotating core can lead to a corresponding periodicity only if the core is 
asymmetric. Hence it seems plausible to attribute periodicities in neutrino, irra- 
diance, and decay-rate data in the range 11-12 year -1 to processes involving 
the rotation of a solar core that is not cylindrically symmetric. This raises 
the possibility that the core may also exhibit north-south asymmetry. In this 
connection, it is worth noting that, since the Sun's axis is inclined with respect to 
the normal to the ecliptic, a north-south asymmetry might contribute to the an- 
nual variation found in decay rates ( | Jenkins et ai, 200 9b Fischbach et ai, 2009; 
Javorsc k II et al, 2010[ ). This effect, in association with an annual variation due 
to the eccentricity of the Earth's orbit, may explain why the phases of the annual 
variations in decay data do not agree with what would be expected on the basis 
of a purely orbital effect. This issue will be discussed in detail in a subsequent 
article. 

These results raise questions concerning not only solar physics but also par- 
ticle physics. They also point to the need for new experiments that monitor the 
decay rates of a variety of nuclei with high cadence and high accuracy. 
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